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ABSTRACT 

General relativistic numerical simulations of magnetized accretion flows around 
black holes show a disordered electromagnetic structure in the disk and corona and 
a highly relativistic, Poynting-dominated funnel jet in the polar regions. The polar 
jet is nearly consistent with the stationary paraboloidal Blandford-Znajek model of 
an organized field threading the polar regions of a rotating black hole. How can a 
disordered accretion disk and corona lead to an ordered jet? We show that the polar jet 
is associated with a strikingly simple angular-integrated toroidal current distribution 
dirp/dr oc r~^/*, where I^{r) is the toroidal current enclosed inside radius r. We 
demonstrate that the poloidal magnetic field in the simulated jet agrees well with the 
force-free field solution for a non-rotating thin disk with an r~^/^ toroidal current, 
suggesting rotation leads to negligible self-coUimation. We find that the polar field 
is confined/coUimated by the corona. We also study the properties of the bulk of 
the simulated disk, which contains a turbulent magnetic field locked to the disk's 
Keplerian rotation except for rapidly rotating prograde black holes {a/M > 0.4) for 
which within r < 3GM/c^ the field locks to roughly half the black hole spin frequency. 
The electromagnetic field in the disk also scales as r"^/'', which is consistent with some 
Newtonian accretion models that assume rough equipartition between magnetic and 
gas pressure. However, the agreement is accidental since toward the black hole the 
magnetic pressure increases faster than the gas pressure. This field dominance near 
the black hole is associated with magnetic stresses that imply a large effective viscosity 
parameter a ~ 1, whereas the typically assumed value of a ^ 0.1 holds far from the 
black hole. 

Key words: accretion disks, black hole physics, galaxies: jets, gamma rays: bursts. 
X-rays : bursts 



1 INTRODUCTION 

Black hole accretion is one of the most powerful sources of 
energy in the universe. A substantial fraction of the gravita- 
tional binding energy of the accreting gas is released within 
tens of gravitational radii from the black hole, and this en- 
ergy supplies the power for a variety of astrophysical systems 
including active galactic nuclei, X-ray binaries, and gamma- 
ray bursts. Elucidating the processes that take place in the 
central regions of black hole disks is obviously crucial if we 
wish to understand the physics of these energetic objects. 

Magnetized, differentially-rotating accretion disks ex- 
hibit the magneto-rotational instability (MRI) and magne - 
tohydrodynamic turbulence (|Balbus fc Hawlevlll99ll . [l998l ). 
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which generate large spatio-temporal variations in all fluid 
quantities and strong correlations between fluid quan- 
tities. Recent general relativistic magnetohydrodynamic 
(GRMHD) simulations of black hole accretion systems 
have begun to resolve these processes and have revealed 
a flow structure that can be decomposed into a disk, 
corona, disk wind, a nd highly magnetized polar region 
that contains a jet (IDe Villiers. Hawlev. fc KrolikI l2003l : 
iMcKinnev fc Gammiell2004l ). As expected, the simulations 
show complex time-dependent behavior in the disk, corona, 
and wind. Surprisingly, however, the polar regions of the flow 
are found to have a simple structure with a nearly force-free, 
time-steady Poynting-dominated jet |McKinnev fc Gammid 
l2004l : iHawlev fc Krohkll2006l : lMcKinney||2006d) The numer- 
ical solution here is quantitatively co nsistent with the rel- 
ativistic force-free m odel proposed by 'Blandfor d fc Znaiekl 
(il97S), hereafter BZ (|McKinnev fc Gammia,20o3 r 
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The primary question this paper explores is the follow- 
ing: how can a turbulent accretion disk lead to a nearly 
stationary, collimated and ordered Poynting-dominated jet? 
In simple force-free models of the disk-jet coupling, like the 
one developed by BZ, one finds stationary solutions for a 
fixe d toroidal current and angu lar velocity of the disk (see, 
e.g. iBlandford fc Znaiekl [iQTtI ). Since in GRMHD simula- 
tions the poloidal field threading the black hole is simple 
and nearly stationary, this implies that the toroidal current 
must also be simple and nearly stationary. However, it has 
not yet been known or understood what radial dependence 
would be chosen by turbulent accretion flows driven by the 
MRI. In order to place the GRMHD simulation results in 
the context of analytical models of jets and winds that treat 
the disk as a equatorial boundary condition, our first ob- 
jective is to determine the radial dependence of the toroidal 
current, angular velocity of the plasma, and angular velocity 
of the magnetic field. 

Given the radial distribution of the toroidal current and 
angular velocities of the plasma and field, one can generate 
force-free models of the jet that approximate the disk as an 
infinitely thin rotating conductor. General relativistic force- 
free solutions of this kind are obtained and then compared 
to GR MHD simulations and the mo dels of BZ in a foUowup 
paper l|McKinnev fc Naravanlbood ). 

Our second objective in this paper is to determine the 
angular distribution of toroidal currents and the flow pattern 
of poloidal currents in the accretion flow. The location of the 
toroidal currents helps identify the role played by the weakly 
magnetized corona in confining the highly magnetized jet. 
For example, the corona provides forces that balance the 
forces due to the strong poloidal field gradients (associated 
with strong toroidal currents) at the boundary between the 
magnetized jet and corona. Hence, the corona can be un- 
derstood as required to confine the magnetized jet. Also, a 
stationary model must have poloidal currents that close like 
a circuit, yet it has not been known where such poloidal 
currents fiow in MHD turbulent disks around black holes. 
We study the distribution of poloidal currents that are as- 
sociated with the outgoing power of the jet/wind in order 
to establish where the poloidal currents flow and to resolve 
the issue of current closure. 

Our third objective is to check if either of the two 
magnetic field geometries described by BZ, viz., the split- 
monopole and the para boloi dal geometries (originally dis- 
covered bv Miche]| |l973l and lBlandfordlll976l . respectively), 
is a good description of the jets found in our GRMHD sim- 
ulations. We find that neither model is satisfactory. Instead 
we identify a third model, in between the other two and 
close to the paraboloidal model, which agrees surprisingly 
well with the simulations as long as the jet is nearly force- 
free. 

The toroidal current in the split-monopole and 
paraboloidal solutions scale with radius as dl^/dr oc 
r~ , r~ , respectively, whereas the toroidal current in 
the GRMHD simulations is found to scale as r"^''*. In- 
terestingl y, the latter scaling is i dentical to that pro- 
posed bv IBlandford fc Pavnd l|l982l ) (BP), who developed 
a non-relativistic self-similar magnetohydrodynamic (MHD) 
model of disk winds by assuming that the sound speed and 
Alfven speed in the disk scale similarly with radius. Our 
fourth objective in this paper is to establish how the disk 



magnetic field strengths, sound speed, Alfven speed, plasma 
speed depend on radius in order to determine whether the 
agreement between the GRMHD simulations and the BP 
model has a deep physical significance or is merely a co- 
incidence. The answer appears to be the latter in the sense 
that the assumptions made by BP are broken near the black 
hole. More importantly, the field threading the disk is dis- 
organized and the disk wind is thermally-driven instead of 
behaving like a "bead on a wire" as in the BP model. 

Prior studies of the BZ power output suggested that 
the black hole power output should be too small compared 
to the disk power output and too small to account for the 
most powerful radio sources l|Ghosh fc Abramowic3 Il997l : 
iLivio. Ogilvie. fc Pringlelll999l '). Such studies assumed that 
the effective viscosity parameter a was as determined in non- 
relativistic simulations and assumed the field strength near 
the black hole was set by sub-equipartition arguments. Our 
flnal objective is to determine the magnetic a viscosity pa- 
rameter as a function of radius within the disk. Both of their 
assumptions end up not applying near the black hole. 



PAPER OUTLINE 

In section [2l we discuss the origin of the ordered poloidal 
fleld in GRMHD simulations. We show that the angular- 
integrated toroidal currents in the turbulent accretion disk 
follow a simple power-law behavior. We discuss the angular 
structure of the toroidal currents and the flow of poloidal 
currents in the accretion flow. We discuss the field angular 
velocity in the transition region between the accretion disk 
and black hole. In section |3l we study the GRMHD accre- 
tion flow in order to extract other electromagnetic proper- 
ties, such as the magnetic fleld strength as a function of 
radius. We test the assumptions of BP against our GRMHD 
numerical models and study the electromagnetic stress that 
leads to an enhanced angular momentum transport near the 
black hole. In section |4l we discuss the limitations of our 
calculations. Finally, in section [S] we discuss our results and 
conclude. 

In appendix [X] we summarize the GRMHD equations 
of motion and point out the reduction to the force-free set 
of equations. In appendix [B] we show how to obtain force- 
free solutions in Schwarzschild and flat spacetimes for an 
arbitrary current sheet at the equator, and we discuss how 
the disk currents are integrated to obtain a toroidal current 
density as a function of radius. We also give three example 
solutions corresponding to the split-monopole, paraboloidal, 
and our new self-similar solution. 



UNITS AND NOTATION 

The units in this paper have GM = c = 1, which sets the 
scale of length {r^ = GMji?) and time { tg = GMje" ). The 
horizon is located at r = r+ = rg(l + sj\ — (a/M)^)). For 
a black hole with angular momentum J = aGM/c, a/M is 
the dimensionless Kerr parameter with —1 ^ a/M ^ 1. In 
order to obtain a density for a given mass accretion rate, one 
requires the fleld as a function of bl ack hole sp i n given by 
GRMHD models such as described in lMcKinnevI l|2005al lbll3. 
l2006d ). The mass scale is determined by setting the observed 
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(model-dependent measured or inferred) mass accretion rate 
(Mo) equal to the accretion rate through the black hole hori- 
zon as measured in a simulation. So the mass scale is set 
by the mass accretion rate (Mo) at the horizon, such that 
Mo[r = r+Jtg/rg and the mass scale is then just 



PO,disk = 



,fcr^ = Mo[r = r+]tg. 



The results of the simulations can be applied to 
any astrophysical system once the value of po,disk is es- 
timated. For example, a coUapsar model with M = 
O.IMqs"^ and M » 3Mq ha s po,disk ~ 3.4 x lO^^gcm"^ 
(|MacFadven fc WooslevI Il999l ). M87 has a mass accretion 
rate of M ~ lO '^Mp) yr~ ^ and a black hole m ass of 
M « 3 X IO^Mq (|HoI |1999I : iRevnolds etahl ll996D giving 
po,disk ~ 10~^®gcm~^. GRS 1915-1-1 05 has a mass accre- 
tion r ate of M ~ 7 X 10~'^M(7) yr ~ ^ (iMirabel fc Rodrigued 
1 1994 iMirabel fc Rodn'guej 1 19991: iFender fc Bellonil |2004| ) 
with a mass of M ~ 14M0 ( Greiner et al.l l20o" ) (but see 



iKaiser et aLll2004h . This gives pp. disk ~ 3 x IQ -'^gcm"^. 

The notation follows Misner et al.l l| 19731 ) and the sig- 
nature of the metric is h 4-4-. Tensor components are 

given in a coordinate basis. The components of the ten- 
sors of interest are given by g^v for the metric, F'^" for 
the Faraday tensor, F for the dual of the Faraday, and 
T'^" for the stress-energy tensor. The determinant of the 
metric is given by ^/—g = Det{g^,^). The field angular fre- 
quency is Of = Ftr/Fr4, — Fte/Fe^. The magnetic field 

* it 

can be written as B = F . The poloidal magnetospheric 
structure is defined by the (/!>-component of the vector po- 
tential Aijy. A stationary, axisymmetric current system is 
defined by the current density J and the enclosed (from 
the pole to some point) poloidal current [B^ = F^t)- The 
elec tromagnetic luminosity i s L = — 27r d6T^^ '^\ r^ sin 6. 



See iGaminie et all ll2003ah: iMcKinnev fc Gammig l|2004h : 
iMcKinnevI (|2004l . l2005bl ll l2006al l for details on this stan- 
dard notation. 



2 THE ORGANIZED POLAR FIELD 

In this section we discuss the origin and nature of the orga- 
nized field threading the black hole. We first review some rel- 
evant results from GRMHD simulations of accretion flows. 
We then demonstrate that at large radii the jet from the 
black hole is electromagnetically pure, while the disk wind 
is dirty. Next, we show that the angular-integrated toroidal 
current over the accretion flow is simple and has a power- 
law radial dependence, which is associated with the simple 
organized polar field. We analyze the angular structure of 
these toroidal currents to locate the toroidal currents (or 
equally the large poloidal field gradients) . We show how the 
jet and disk wind are associated with poloidal currents that 
close in an electric circuit. Next, we determine the radial 
dependence of the disk-averaged plasma and field angular 
velocities, which have a simple behavior. This also eluci- 
dates how strongly the field in the disk couples to the rota- 
tion of the black hole. Finally, we compare the polar field in 
the GRMHD simulations to the non-rotating force-free field 
solution that emerges from the same power-law toroidal cur- 
rent as in the GRMHD simulations. The agreement found 
between these models demonstrates that the polar jet is ac- 
curately described by a force-free model and that rotation 
plays a negligible role in self-coUimating the polar jet. Ro- 



tation may still play a role in indirectly coUimating the jet 
via the rotation of the surrounding coronal wind. 



2.1 Review of Prior GRMHD Simulation Results 

In this section, we review some relevant results from 
GRMHD simulations, w here we start the di scussion with 
the fiducial model of McKinnev fc Gammij f2004). Fig- 
ure [T] shows the time-averaged field geometry for this fidu- 
cial model. The black hole has a spin of a/M — 0.9375, 
which is close to the equilibrium va lue of a/M ~ 0.92 
(jGammie. Shapiro, fc McKinnevll2004l ). 

The simulated accretion flow is dominated by a tur- 
bulent hydromagnetic dynamo driven by the MRI, which 
drives an increase in poloidal field strength whose saturated 
magnitude is insensitive to the initia l poloidal field strength 
toe Villiers. Hawle v. fc KrolikI l2003l : iMcKinnev fc Gamniiel 
[2004). The dyna mo eventually d ies out unless modelled in 
three dimensi ons ([Cowlinel 1934t). For the two-dimensional 
simulations of IMcKinnev fc Gammij |2003), measurements 
are made only during the turbulent period and their results 
are consistent with three-dimensional simulations. 

The polar region contains a well-ordered field whereas 
t he disk appears to have a t urbulent and diso r dered field 
dMcKinnev fc Gammij|2004l : iHirose et aLll2004 iMcKinnevI 
l2005a|). Between these two regions is the corona which, in 
a time-averaged sense, contains only weak disordered fields. 
In the figure, the field that appears to come from the disk 
does not reach large distances, where as the organized fi eld 
in the polar region reaches large radii (|McKinnevll20063 ) . 

It is important to note that the organized polar field ge- 
ometry shown in Figure [1] is generic for simulations that are 
initialized with a well-organized poloidal field with or with- 
out net flux. Here an "organized" poloidal field means a field 
geometry with few poloidal polarity changes, but it could 
be initially contained entirely within the disk. The quasi- 
stationary structure of the accretion flow and flux threading 
the black hole or disk otherwis e depends little on the initia l 
fleld strength or geometry Q (|McKinnev fc Gammid |2004| : 
iHirose et al]|2004l '). 

Cases when the black hole does not end up with an or- 
ganized field include an initially purely toroidal field (no 
toroidal currents so no poloidal field) (|De Villiers et all 
l2005al ) or highly tangled field (McKinney fc Gammie, in 
prep.). With a purely toroidal field, no jet is produced. With 
a highly tangled field, the Poynting-dominated jet does not 
form since it is continuously contaminated with disk ma- 
terial, and so the corona and coronal wind dominate the 
entire polar region. Even when the disk is initially threaded 
by net fiux, the quasi-stationary corona still only contains 
weak disordered field. 

The lack of organized field in the corona or threading 



The codes used in the studies mentioned above preserve the 
solenoidal constraint (V ■ B = 0) to machine accuracy, and so the 
codes preserve the net poloidal flux in the system. Hence, the net 
radial flux remains constant and the 6 flux only changes by losing 
flux into the black hole or through the outer radial boundary. 
Codes that violate this property might generate large artificial 
net flux and then regions of organized flux could not be trusted. 
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the disk might be understandable as a result of the inabil- 
ity of flux to be simply advected radially inward if the ra- 
tio of viscosity to magnetic diffusivity (magnetic Prandtl 
number Pm) is such that P,n < H/R, where H/R is the 
pressure sca le-height per unit rad ius of the disk (see discus- 
sion in, e.g.. lLivio. Qgilvie. fc Prin gle 1999. sections 3.1 and 
3.2). In thin turbulent-driven accretion flows, the magnetic 
Prandtl number is order unity, and so one does not expect 
to be able to simply accrete net flux from large radii. On 
the other hand, thick advection-dominated accretion flows 
are expected to be able to advect a nonnegligible flux. In 
three-dimensional models, the advection of flux may occur 
in isolated flux tubes an d still build near the black hole 
l|Spruit fc Uzdenskvllioosi ). The physics of how large-scale 
flux can be accreted through disks remains an open is sue 
(see, e.g.. lRevnolds et al.|[200^ : IContopoulos et al.ll2006l ). 

So in the GRMHD simulations, how does the black hole 
become threaded with flux and how is the magnitude of 
the field maintained? Does large-scale fiux simply advect 
inward despite the above-mentioned issues? In the simula- 
tions, the initial large-scale flux around the black hole is 
created as a result of the accretion of equatorial fleld loops 
that do not thread vertically through the disk, thus bypass- 
ing the problem of advecting large-scale flux threading the 
entire disk. The long-term flux is maintained by the 2D ax- 
isymmetric dynamo that drives single polarity poloidal loops 
(corresponding to single polarity toroidal currents) to twist 
in an axisymmetric sense, break through reconnection, and 
interchange around each other within the disk allowing the 
black hole to accrete an arbitrary poloidal polaritj0. Near 
the black hole the electromagnetic field dominates and this 
allows the magnetic flux threading the horizon to grow by 
the attraction of similarly-signed toroidal currents and the 
repulsion of oppositely-signed toroidal currents - a classical 
electrodynamical phenomena. 

The magnetic field strength near the black hole sat- 
urates when material forces of the disk-|-corona can just 
support the magnetic pressure of the polar field. Thus 
one expects the magnitude of the magnetic pressure to 
be in near equipartition with the magnitude of the gas 
pressure fairly clo se to the horizon, which is the case 
l|McKinnev fc Gamm ic 2004). This balance is qualitatively 
similar to the balance that is associated with the "magne- 
tospheric radius" in accreting neutron star systems. 

Also important is the accretion of large-scale polarity 
changes that can overwhelm any prior fleld built-up around 
the black hole. The relevance of accreting large-scale oppo- 
site polarities probably depends o n the astrophysical sys- 
tem (see, e.g.. iNaravan et al.ll2003l) . Reconnection tends to 
convert the pre-existing flux into thermal energy and erase 
the organization of the polar field. When the reconnection 
time-scale is long, the fleld energy associated with the flux 
is difficult to remove once the fiux is in place. 

The magnetized jet is produced as the black hole spin 
and disk rotation create a large toroidal field whose gradient 
launches a significant fraction of accreted fiux out along the 
poles. This process effectively leads to net fiux threading 



^ In a 3D dynamo non-axisymmetric modes would more vigor- 
ously generate such varying polarities and allow direct radial in- 
terchange not allowed in axisymmetry. 




Figure 1. Time-averaged poloidal magnetic field (solid black 
lines) for the fid ucial GRMHD numerical mo del with a/M = 
0.9375 studied bv lMcKinney fc Gammij ||2004| ). Color shows the 
logarithm of the time-averaged rest-mass density (log po ) from 
highest (red) to lowest (blue) densities. The black hole is located 
at (R, z) = (0, 0). Notice (i) the ordered, magnetized, low-density 
jet near the polar axis, (ii) the turbulent, disordered disk in the 
equatorial region, and (iii) the coronal region in between with 
only a weak disordered field. 

the black hole since the flux that reaches larg e distances 
becomes causally disconnected from the disk /Mc KinnevI 
[2006c,). The stronger the black hole spin, the stronger is this 
effect. In a stationary state, the organized polar field threads 
the black hole horizon, but not the inner disk (|Hirose et al.l 
I2OO4I : lMcKinnevl[2005al ) . 

In summary, the black hole naturally becomes threaded 
by organized field when the disk contains a field with few 
large-scale poloidal polarity changes. Early in the simula- 
tion, the field threading the black hole is advected through 
the equatorial region rather than being advected as large- 
scale flux threading the entire disk. The magnetic fleld 
strength near the black hole saturates when a balance 
is reached between the polar magnetic pressure and the 
disk+corona pressure. The disk and black hole rotations 
drive the flux near the black hole to large radii, and this 
effectively leads to a net flux threading the horizon. In a 
quasi-stationary state, large-scale field threads the hole but 
not the inner disk. 



2.2 Power Output of Black Hole and Disk 

In principle, a significant fraction of the electromagnetic 
power may come from a disk wind (iBlandford fc Znaiekl 
and may even dominate the black hole power out- 
put ('Gho sh fc Abramowicj 1 19971 : iLivio. Ogilvie. fc Pringlj 
Ji)99). However, the lack of an organized fleld threading the 
disk suggests that the electromagnetic output of the disk 
may be seriously compromised compared to those simple 
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Figure 2. The angular density of tlie electromagnetic power 
[dP/dd) vs. 9 in radians at four selected radii, r = 
{r+, nsco: lO^'g: 40rg}, fo r the fiducial GRMHD model with 
a/M = 0.9375 described in iMcKinnev fc Gammi3 ||2004|') . The 
thick solid line corresponds to r = 40rg, the outer radius of the 
computational domain for this simulation. Only the black hole 
electromagnetic power output survives at this distance, whereas 
the disk electromagnetic power output has been efficiently con- 
verted into other forms. 



models that have treated the disk as just a boundary con- 
dition. To test these ideas, we compare the electromagnetic 
power output of the black hole and the disk in the fiducial 
[a/M = 0.9375) model shown in Figure [T] 

Let us define the time-averaged angular density of the 
electromagnetic power output as 



2nr 



-T 



(Euy 



(1) 



where _y(^*f)^ ig the radial electromagnetic fiux written 
in a suitable coordinate basis (either Boyer-Lindquist or 
Kerr-Schild coordinates). The time-averaging is performed 
over approximately 8 orbital periods at r = IQvg once the 
fiow reaches a quasi-stationary, turbulent state. Figure [2] 
shows {dP/dO) as a function of 9 at four radii: the hori- 
zon at r = r+, the inner- most stable circular orbit (ISCO) 



at r 



nsco, r 



AQva- Notice how smooth 



and well-behaved the power is near the poles and how erratic 
and disordered it is away from the poles. 

The total electromagnetic power at any radius r is given 

by 



(2) 



Table [T] lists values of P at the same four radii shown in Fig- 
urelH For comparison with th e force-free models discussed in 
iMcKinnev fc NaravanI (|2006| ). the power is normalized such 
that 



P(r) 



P(r) 



((B'-)'|e=o)47rr2c' 



(3) 



where the field is in Gaussian units and is measured on the 
horizon at the poles at the final time of the simulation. 

At the horizon the electromagnetic power is that from 
the black hole with most of the net electromagnetic power 
coming from the polar regions. With increasing radius, elec- 
tromagnetic power from the disk is added. However, this 
power is steadily converted into matter energy so that, by 
r = 40rg, most of it is no longer in electromagnetic form. 
This explains why the power at r = 40rg in Table [T] is not 
much larger than the power at r = r+. 

Table [1] also gives the Lorentz factor F of the jet far 
from the black hole and the half-opening angle 6j of the jet 
(defined as the angle at which (dP /dO) is maximum). These 
results are from iMcKinnevI (|2006cf ) . 

In summary, we find that the electromagnetic power 
at the poles of the black hole remains undiminished out to 
the largest radius shown, whereas the electromagnetic power 
from the disk region decreases outward as it is efficiently con- 
verted into kinetic and thermal energy of the plasma. Sim- 
ple estimates of the power output of the black hole and disk 
did not consider the conversion of electromagnetic power 
into t hermal and material power in the corona and disk 
wind (|Ghosh fc Abramowicj|l997l : lLivio. Ogilvie. fc Pringlel 
1 999j ), and so they may have seriously overestimated the 
power from the disk. Thus, the black hole may dominate 
the electromagnetic power output of accretion systems. Also, 
when the disk is turbulent, the black hole remains the only 
possible clean {h^ /{poc? 
power. 



3> 1) source of electromagnetic 



2.3 A Power-Law Radial Distribution of Toroidal 
Current 

For stationary fiows, the toroidal current directly leads to 
the poloidal field structure. Since the Poynting-dominated 
jet has a simple, stationary poloidal structure, there must 
be simple toroidal currents to support the field. However, 
given the turbulence in the disk and the efficient conversion 
of electromagnetic power from the disk into material power, 
one might assume that the electromagnetic properties of the 
disk would be complex and hard to relate in any simple way 
to the organized polar fiux. This is certainly suggested by 
Figures [T] and [5] 

As discussed in the introduction, the motivation for 
studying the toroidal current comes from simplified force- 
free models that include the accretion disk as an equa- 
torial boundary condition, such as the paraboloidal BZ 
model (|Blandford fc Znaiekl [l973)- In general, models of 
winds and jets typically specify the toroidal current in 
the disk (or stellar surface) and find the corre s pond- 
ing solution (se e e.g. iMichej Il973l: lOkamotd Il974 ll978l: 



BlandfordI Il976l: iBlandford fc Pavii^ Il982l: [Lovelace et al 



19861 : iHevvaerts fc Normaijl 19891: iNitta et al.lll99ll:lLi et al 



199?; Appl fc Camcnzind' T992', T993; 'Beskin fc Pariev 
1993; Contopoulos 1994; Contopoulos fc Lovelace 1994|: 



Contopouloslll995al lbl). We stress that for a stationary solu- 
tion the toroidal current and poloidal field are just different 
languages for the same physics, but our motivation is to see 
if a turbulent GRMHD disk can be related to the vast array 
of models that describe the disk as simply a current sheet. 

The relevant questions are: 1) Are such simple models 
applicable to thick, turbulent, magnetized accretion flows? ; 
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GRMHD Model Electromagnetic Power 



a P[r+] P[nsco] P[r = 10rg] P[r = AOrg] r[r = 5 X lO'^T-g] 6^ [r = 5 X IQ-^rg] 

0.9375 0.324 0.407 0.595 0.374 10 5° 

Table 1. Electromagnetic power (per unit (B^)^ on horizon at poles), Lorentz factor, and half-opening angle. 



2) What radial dependence of the toroidal current is the "cor- 
rect" choice? ; and 3) Does the angular-integrated toroidal 
current predict the shape of the organized polar field? 

Let us consider the angular-integrated (over all angles) 
toroidal current that for a stationary flow corresponds to a 
line integral of the magnetic field around a closed poloidal 
loop. We integrate over all angles to capture currents that 
sometimes rise into the corona and to capture variations 
in the disk-|-corona thickness in time and as a function of 
radius. We consider the angular distribution of toroidal cur- 
rents in the next section. 

The current density (as given by equation (|A11[| ) can 
be integrated to determine the net toroidal current enclosed 
within the volume between ro and r. This invariant current 



IS 



= r r v^(j^dr'de') 

Jrn Je=0 ^ ' 



ro J 9=0 
dh. 



(4) 



where dE^ = e^.^apt" dr°'de'^ , t'' = {1,0,0,0} is the time- 
like Killing vector, and dl /dr is the toroidal current per unit 
radius. Note that the magnitude of the enclosed current is 
set by the value of I^{ro), which is an arbitrary constant 
and set to be the enclosed current in the numerical grid of 
size dr at tq. Only the magnitude of the current density has 
physical significance. Notice that 



dlcf, 
dr 



-g,rde, 



(5) 



where ^—g « r sin 9 far from the black hole. 

In the following we are particularly interested in models 
with a power- law dependence of the current density, i.e.. 



dl^ 
dr 



(6) 



We are motivated by the fact that the split-monopole and 
paraboloidal force-free models both have currents of this 
form, with = and 1, respectively. Example solutions 
are given in appendix iBl 

The solid line in Figure [3] shows the enclosed toroidal 
current /0(r) at a very earl y time (t = 5Qtg) of t h e fidu - 
cial simulation described in iMcKinnev fc Gammid (|2004 ). 
At this time, the system has hardly deviated from the ini- 
tial conditions and there is very little accretion taking place. 
Correspondingly, the current has a fairly complicated depen- 
dence, which primarily reflects the particular initial condi- 
tions chosen for this simulation. 

Figure 3] shows the enclosed toroidal current at a later 
time {t — lOOOtg) when the accretion flow is highly turbu- 
lent and has reached a quasi-steady state. We see that the 




Figure 3. Vertically-integrated enclosed toroidal current as a 
function of radius. Solid line shows the current in the GRMHD 
numerical model with a/M = 0.9375 at t = 50tg (nearly the 
initial conditions). The dotted line is the u = monopole solu- 
tion, the long-dashed line is the u = I paraboloidal solution, and 
the short-dashed line is the !^ = 3/4 solution. Clearly the initial 
conditions do not follow any simple power-law behavior. 



current distribution has changed dramatically from its ini- 
tial distribution. More importantly, the current profile looks 
smooth and simple. Figure [5] shows the enclosed current as 
time-averaged over the period t = 500tg to 1500tg (roughly 
a turbulent time scale at r = 40rg). The result is very simi- 
lar to that shown in Figure |4l except that the current looks 
even smoother because of the averaging. 

In Figures 1 3151 we show for comparison the enclosed 
currents corresponding to the split-monopole {u = 0, dot- 
ted line) and paraboloidal {v — 1, long-dashed line) mod- 
els. It is clear that neither of these power-law models is a 
good representation of the enclosed current in the steady 
state GRMHD model. On the other hand, the short-dashed 
lines, which correspond to a power-law model with v = 3/4, 
describe the quasi-stationary GRMHD results surprisingly 
well. This particular model is associated with a radial de- 
pendence of the current density of the form dl^/dr oc r"®/*. 

A few interesting conclusions can be reached from these 
results: 1) The GRMHD model is nearly coincident with the 
i/ = 3/4 model; 2) The GRMHD model is not consistent 
with the split-monopole or paraboloidal models; 3) Despite 
the complicated nonlinear turbulence, the currents in the 
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Figure 4. Similar to Figure |3] except tliat tlie solid line shows 
the result for the GRMHD numerical model at t = lOOOtg , during 
a period of strong sustained disk turbulence. By the action of 
the magneto-rotational instability, the GRMHD model's toroidal 
current has redistributed itself to closely follow a i/ = 3/4 power- 
law dependence. 




0.01 - 



disk are simple not only on average, but at each moment in 
time. 

The simple = 3/4 current distribution described 
above is entirely consistent with the smooth field distribu- 
tion seen in the evacuated polar region in Figure [T] It is also 
consistent with the fact that the Poynting-dominated jet is 
nearly stationary and nearly resembles BZ's paraboloi dal so- 
lution (jMcKinnev fc Gammidliooj : lMcKinnevir2006cl ) . Fur- 
thermore, the fact that the best fit is obtained for = 3/4 
rather than u = 1 explains why the large-scale field lines in 
the jet are nearly paraboloid al but somewhat less coUi mated 
in the GRMHD simulations (|McKinnevll2005d . 120063 '). Sec- 
tion [2]7] discusses this point further. 

The above results are roughly independent of the ini- 
tial field geometry as s umed in the GRMHD simulations. 
iMcKinnev fc Gammid l|2004l ) considered different initial 
conditions such as multiple magnetic loops in the initial 
torus, a net vertical field, and loops of alternating poloidal 
direction. Once these simulations have reached a quasi- 
steady state, each closely follows a power-law toroidal cur- 
rent density with v — 3/4. This is despite the fact that the 
particular angular structure of the toroidal current varies 
considerably. For example, for thin disks with alternating 
polarity of multiple field loops, there is no strong poloidal 
field in the funnel, but the toroidal current still maintains 
the u = 3/4 power-law dependence. 

We also investigated the dependence of the toroidal 
current distribution on the black hole spin. Once again we 
find that, for aU a/M ranging from -0.999 to +0.999, the 
toroidal current settles down to a = 3/4 power- law dis- 
tribution once the fiow reaches quasi-steady state. This is 
despite significant changes in the enclosed poloidal current 
as shown in Figure [T] since for slowly spinning black holes 
there is negligible poloidal current flowing in the funnel re- 
gion. 

How robust are these results to changes in the mass 
distribution and disk thickness? The flducial model studied 
has an initial mass distribution of a hydrostatic equilibrium 
torus with a constant speciflc angular momentum such that 
H/R ~ 0.3. An alternative initial mass distribution we tried 
has a quasi-equilibrium Keplerian disk with a Gaussian ver- 
tical mass distribution with H/R ~ 0.3. Once the turbulence 
reaches the nonlinear phase, this alternative model also has 
a toroidal current closely matching the i/ = 3/4 profile. 

We have also studied a thin (i.e. small H/R) Keplerian 
disk with an initial Gaussian vertical distribution with an ad 
hoc cooling model to keep H/R ~ 0.05, and we find that this 
model also obeys the = 3/4 toroidal current distribution. 
The only qualitative change is that there is a larger variance 
around the v = 3/4 solution. However, the solution is still 
quite different from the monopole or paraboloidal solutions. 



Figure 5. Similar to Figure [J] except the enclosed toroidal cur- 
rent has been time-averaged over several turbulent dynamical 
times at r ~ lOrg. Evidently, the simple u = 3/4 power-law be- 
havior for the toroidal current holds for both the time-averaged 
quasi-stationary turbulent state and at every moment in time 
once the turbulence has grown to saturation. 



2.4 Angular Structure of Toroidal Current 

Of course, the angular-integrated toroidal current does not 
reveal the angular location of the currents. By integrating 
equation Q over a finite radial range for each 6, one can 
identify structures in the accretion flow with the toroidal 
currents. While the instantaneous value of I(i,{9) is very 
oscillatory and difficult to interpret, the running integral 
ilo Iti>{9')d6') is relatively smooth so we focus on this quan- 
tity. 
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We consider two radial sections of the accretion flow. 
One radial range considered is r = r+ to r = 3rg, while 
the other radial range considered is r = r+ to r = lOrg . Fig- 
ure|6]shows the radially integrated, a ngular running integral 
of Id>( d) for the fiducial simulation in lMcKinnev fc Gammig 
((2OOJ). For this GRMHD simulation, the funnel jet extends 
from the poles inward by about 1/2 to 1 radian depending 
upon the radius. Just beyond this range at the corona-funnel 
interface, the toroidal current changes sign and nearly an- 
nihilates itself in an integral sense. This toroidal current is 
associated with the polar field. In the corona there is no or- 
ganized field and the magnetic field is in equ ipartition with 
the gas pressure (jMcKinnev fc Gammiell2004l '). We find that 
across this coronal-funnel interface there is force balance be- 
tween the organized field in the polar region and the equipar- 
tition, disorganized, weakly magnetized, hot plasma in the 
coronal region. If self-coUimation is weak within the polar 
jet, then the corona is responsible for confiriirig/collimating 
the polar jet. The importance of rotat ion of the polar field is 
tested in later sections and in iMcKinnev fc NaravanI (|200^ . 

The next jump in the toroidal current is at the corona- 
disk interface. Within the disk the toroidal current oscillates 
with a slightly non-zero toroidal current. For the integrated 
radial range of r = r+ to r = Wrg, the angular integrated 
distribution of the toroidal current is roughly given by 



I{e')de' oc (l-cos(6i)). 



(7) 



although there is significant sub-structure. This fit is shown 
in Figure [6] 

While changes in the initial field geometry do not affect 
the radial distribution of toroidal currents with dJ^/dr oc 
r~^^^, the angular location of the toroidal currents (and so 
the location of poloidal field gradients) depends sensitively 
on the initial field geometry. Simulations that start with a 
highly disorganized field with many loop of different polarity 
lead to no organized flux in the poles and then the toroidal 
currents are primarily located inside the disk-corona inter- 
face. 



2.5 Poloidal Currents and Current Closure 

For stationary flows, the toroidal field is equally described 
by the poloidal current, and these are related to the power 
output of the jet. At large radii, the radial electromagnetic 
power output is given by -B oc ilpB^B^, where Of is the 
field rotation frequency, is the radial field strength, and 
B^ is the "polar enclosed poloidal current" {B^ « RB"^). 
In a stationary, axisymmetric fiow, Q.f and B^ are constant 
along field lines labelled by the vector potential (also 
referred to as the fiux function ^ or stream function). Thus, 
the poloidal current is an indicator of the electromagnetic 
power output per unit poloidal field strength, and the closure 
of this poloidal current is best considered in a plot of B^ vs. 

Figure [7] shows the time-averaged value of B^ vs. the 
time-averaged value of A^ at r = lOrg. The values within 
A^ < 0.2 correspond to the region inside the funnel that 
contains the Foynting-dominated jet. The value of A^ « 0.2 
corresponds to the transition between the force-free fun- 
nel and the corona. The value of A,f, « 0.21 is at the 
transition between the corona and disk, where the increase 




1 2 

0[radians] 

Figure 6. Both panels show the running integral of the toroidal 
current I^. Top panel corresponds to the integral over the radial 
range from r = r+ to r = 3rg. Bottom panel corresponds to 
the integral over the radial range from r_|_ to r = lOrg, where a 
rough fit is shown by the dashed line. The vertical dotted lines 
show the location where b'^/{poc^) = 1 (i.e. the funnel-coronal 
interface). The overall toroidal current is seen to be distributed 
over a scale-height of the disk, but the funnel-corona and corona- 
disk interfaces harbor significant toroidal currents. The polar field 
is associated with the toroidal current within the corona near the 
funnel. 



in the poloidal current is largest. Within the disk the en- 
closed current remains constant until the equator is reached 
and there is a strong return current at « 0.26. No- 
tice that both hemispheres have been shown to demon- 
strate that despite the time-averaging, the equatorial region 
is still asymmetric while the jet region remains symmetric. 
Such a plot can be used to compar e GRMHD accretion sys- 
tems to pulsar sy stems (see, e.g.. IContopoulos et al.|[l999l : 



sys 

iMcKinnevI 2006b[) and simila r black hole force-free systems 
( McKinnev fc NaravanI I2OO6I I. Despite the turbulence, the 
dependence of B^ on A^ remain surprisingly simple. Simu- 
lations that start with a highly disorganized field lead to a 
more complicated poloidal electric circuit that is distributed 
over all angles. 



2.6 Turbulence Leads to Simple Field Angular 
Velocity 

For stationary flows, the angular velocity of the field lines 
(Of) is an important quantity that determines the toroidal 
field geometry and determines the electromagnetic power 
output of a jet or wind. Also, a comparison of Hp and the 
fluid angular velocity Q can reveal how well-coupled the mat- 
ter is to the field. In principle, large scale fields can develop 
in the disk and the plasma might slip arbitrarily along the 
field lines, leading to $7^ 7^ O. However, if there is strong 
turbulence, it would lead to a significant random component 
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Figure 7. Poloidal current enclosed from the pole (B^) to the 
point given by the vector potential {A^) for the surface given 
by spherical polar radius of r = lOrg. The value of < 0.2 
is inside the jet. The value of A^ 0.21 is at the disk-corona 
transition. The value of A^j, Ki 0.26 is at the equator. The plot 
shows that the poloidal current increases inside the jet, but is 
more strongly enhanced at the corona-disk interface. The current 
closes by passing through the disk in a somewhat distributed 
current sheet. 
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Figure 8. Angular frequencies of the plasma (< Q >, solid line), 
the field {{Qp), dotted line), the field using a different averag- 
ing procedure (((n^)), short dashed line), and a ZAMO observer 
(^ZAMO I long-dashed line) for a GRMHD accretion disk simula- 
tion with a/M = 0.9375. All frequencies are plotted in units of 
the local Keplerian angular frequency (Ok) and shown in Boyer- 
Lindquist coordinates. Note that < Q >= Ozamo oil the horizon, 
as required in Boyer-Lindquist coordinates. At large radii turbu- 
lence locks the field to the plasma, while at small radii the field 
rotation is locked to the spin of the black hole. 



to the field and the plasma would be unable to slip as much. 
Simple models of the accretion disk assume an arbitrary 
value of Of in the disk, while here we establish which radial 
dependence of Of is motivated by GRMHD simulations of 
turbulent accretion flows near rotating black holes. 

We first investigate the behavior of Of and O as a func- 
tion o f radius in the fiducial model of lMcKinney fc Gammiel 
l|2004l l. One needs to choose some method to space- and 
time-average Of in order to obtain its radial distribution. 
A poor choice would be to directly average Of = Eg / B"^ 
itself because it is highly oscillatory and is undefined at 
positions where the radial component of the field mo- 
mentarily vanishes. Thus, we consider the ratio of space- 
and volume- averaged quantities to obtain a mean angular 
velocity (Of) = (Eg) / (B^ ) . For each radial shell this quan- 
tity is volume- averaged over a disk scale-height and time- 
averaged over approximately 8 orbital periods at r = lOr^ 
once the fiow has reached a quasi-stationary, turbulent state. 
An alternative time-averaging is performed using the abso- 
lute value of each composite quantity {\Ee\ and IB"^]) giving 
{(Of)) = (ISsD/dSn). 

Figure [8] shows the radial profiles of both forms of Of 
normalized by the local Keplerian angular velocity Ok . Both 
methods of averaging Of lead to similar results. The plot 
also shows the angular velocity of the plasma per unit Keple- 
rian (0/Ok) and the angular velocity of the zero angular mo- 
mentum observer (ZAMO) per unit Keplerian (Ozamo/Ok). 
We use Boyer-Lindquist coordinates, with Ok = l/{r'^^^+a). 
Note that Ozamo = Oh ~ a/(2r+) on the horizon. 



As Figure[8]shows, O ~ Of ~ Ok over much of the disk 
for radii r > 2r+. At these radii, there are no large-scale 
fields for the plasma to slip along, so that the plasma and 
field are locked together in a turbulent mixture. However, 
Of becomes somewhat sub- Keplerian near the horizon as 
the field lines become more ordered. Thus, while the plasma 
is forced to corotate with the black hole at the horizon, the 
field lines rotate slowe r with Of > Oh/2. As discussed in 
iMcKinnev fc Gammiel ll2004i'l this behavior for Of is consis- 
tent with the Gammig (199y) inflow model of the plunging 
region. For r > '20rg, the angular velocities deviate from a 
simple behavior because the solution still depends on the 
initial conditions. 

As in the case of the toroidal current distribution, the 
results described here for the angular velocity of the field 
lines and the plasma are quite robust and are independent 
of the assumed initial field geometry, mass distribution, or 
disk thickness. However, the black hole spin has a dramatic 
qualitative effect. 

For spins from a = —0.999 to a = 0.999, all the models 
have the plasma locked to the field at large radii. However, 
close to the black hole there is a qualitative change in the 
results around a/M ~ 0.4. For a/M < 0.4 we find that 
Of ~ Ok ~ Oh even at the horizon because the disk dom- 
inates over the black hole. However, for a/Af > 0.4 we find 
that Of ~ Oh/2 near the horizon as the black hole domi- 
nates over the disk. This transition at a/M ~ 0.4 is consis- 
tent with the fact that there is a qualitative change in the 
energy output of a black hole-disk-jet system at a/M ~ 0.36 
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|Lil l2000l : lMcKinnevll2005al ) , when the Keplerian angular ve- 
locity of gas at the ISCO is equal to the angular velocity 
of the black hole. It is also consistent with the fact that 
there is negligible (or negative) electromagnetic energy ex- 
tracte d from the black hole for spin s a/M < 0.4 for thick 
disks (jMcKinnev fc Gammia.2004 : McKinnevl [2005a'l. Since 
dP/ dd oc Q,f{0,f — 0,h ) on the horizon, we expect dP/ dO ^ Q 
since Q.f ~ when a/M < 0.4. 



2.7 Comparison between GRMHD field and 
= 3/4 force- free field 

We have shown above that the toroidal current and the 
field angular velocity JIf in GRMHD simulations are well- 
behaved and so easy to model. One expects the current dis- 
tribution to be consistent with the nearly force-free funnel 
field geometry if the funnel field is setup directly by those 
currents. Here we test the robustness of the association be- 
tween the power-law index of the angle-integrated toroidal 
current and the polar field geometry. 

We consider a simplified problem in which neither the 
black hole nor the disk rotates. We replace the disk with a 
current sheet with dl^/dr cx r"^''* (corresponding to = 
3/4) at the equatorial plane and assume that the rest of 
the volume is in force-free equilibrium. In Appendix [BJ we 
describe how to obtain force-free solutions in Schwarzschild 
and fiat spacetimes for an arbitrary current sheet at the 
equator with no rotation. 

Figure [9] shows the field geometry from a time- 
dependent GRMHD numerical model overlayed with the 
force- free field corresponding to v = 3/4. We see that there 
is a reasonably good agreement between the two models in 
the funnel region, where the GRMHD solution is Poynting- 
dominated and force-free. The small difi'erences between the 
models could be due to (i) residual weak time-dependence in 
the GRMHD solution, and (ii) additional coUimation in the 
inner regions of the jet due to either rotation (which is ig- 
nored in t he force-free solution plotted here) or coronal pres- 
sure (see iMcKinnev fc NaravanI |2006|) . Also, at very small 
angles the jet is more paraboloidal. However, the differences 
are small. In fact, the agreement between the GRMHD nu- 
merical models and the ;^ = 3/4 force- free model is found 
to be goo d out as fa r as r ~ 10"^rg in the large scale simu- 
lations of lMcKinnevI (|2006cl ). Beyond this radius the inertia 
of the matter becomes nonnegligible in the GRMHD model 
and the force-free approximation is no longer applicable. 

Although the overlay comparison shown in Figure [9] is 
impressive, we caution that the force- free field shown here 
is for the non-rotating case whereas the GRMHD solution 
corresponds to both a spinning black hole and a spinning 
di sk. We discuss the effec t of ro tation on force-free solutions 
in lMcKinnev fc NaravanI (120061 ). 



3 ELECTROMAGNETIC PROPERTIES OF 
THE DISK 

InlMcKinneyl (|2006d ). we presented a detailed study of the 
electromagnetic properties of the Poynting-dominated jet by 
providing both a qualitative description of the jet and many 
quantitative measures of the jet. For example, radial scalings 
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Figure 9. Overlay of the time-averaged poloidal field from the 
GRMHD numerical model (black) and the analytical u = 3/4 
solution described in Appendix B. Only the portion of the v = 3/4 
solution that overlaps the funnel region of the GRMHD model is 
shown. Evidently, a current distribution with u = 3/4 leads to 
a force-free field consistent with the funnel field in the GRMHD 
numerical model. 

within the jet helped isolate the mechanism of jet formation 
and jet stability. 

In this section we first determine the radial scalings of 
the magnetic field components in the disk-|-corona. We con- 
nect the se scaling s to qu asi-analytic models of the accretion 
fiow by iGammij (|l999l ) that describe the disk inside the 
ISCO. Given the simplicity of the time-averaged fields in 
the disk, we suggest that such models might be extendable 
to include the disk beyond the ISCO. Next, we compare the 
GRMHD simulations to simple Newtonian disk models tha t 
have a disk wind (such as that bv lBlandford fc Pavne|[r982l ). 
While such models have a similar magnetic field scaling in 
the disk, the assumptions they make do not hold in the 
GRMHD simulations. Finally, we discuss the role of mag- 
netic stresses near the black hole and compute the effective 
magnetic a viscosity parameter, which near the black hole 
rises to order unity. 

3.1 Radial Dependence of Disk Magnetic Field 

In this section, we focus on the radial dependence of the 
electromagnetic properties in the bulk of the disk+corona. 
Thus, this study excludes the direct properties within the 
jet. However, the mechanisms for disk-jet coupling may be 
better constrained by knowing the radial scalings within the 
disk-|- corona. 

We consider the space-time average of the absolute 
magnitude of the field strengths. The volume-averaged 
value is found for each radial shell over the disk-fcorona, 
which includes only the flow that has b^/po < 1 
(|McKinnev fc Gammiell2004h . Thus the highly magnetized 
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Figure 10. From above, the four solid lines correspond to B'^, the 
comoving field strength |6j, C", and _B® for a GRMHD accretion 
disk simulation with a/M = 0.9375. The magnetic field strengths 
roughly follow power-law behaviors close to the black hole. The 
solution at r > lOrg is still dependent on initial conditions, so 
this region is excluded from the fitting procedure. 



jet is explicitly excluded. For any quantity B within the 
disk+corona we compute 



B = 



IB, 



I' 



J B sin ed6 

J sin ede 



(8) 



which clearly preserves the radial dependence of the quantity 
in question. 

As before, the time-average is computed over the tur- 
bulent period of accretion. The field strengths are given 
in Gaussian units and norm alized by the rest-m ass density 
within the disk as done in iMcKinnevI (|2006d ). such that 
given an estimate of the density of the disk near the black 
hole or the mass accretion rate near the black hole one 
can convert to physical units. Figure [10] shows these field 
strengths as a function of radius; where for r > lOvg the 
initial conditions still contribute to the solution and so the 
field there is not included in the fitting procedure. 

For models with any black hole spin, we find that the ra- 
dial dependence of the comoving field and toroidal strengths 
roughly follow 



0.5 ( — 



(9) 



as found in lMcKinnevI |2005a). Since the toroidal field domi- 
nates within the disk, the value of the lab-frame toroidal field 
(_B'^) is similar in magnitude and follows a similar depen- 
dence except inside the ergosphere within r < 2rg where the 
coordinate efTects of frame-dragging are strong for rapidly 
rotating black hole models. The angular field roughly follows 



IB" 



\/47rpo,diskc2 



0.08 



(10) 



Thus the comoving, toroidal, and angular field nearly follow 
the r~^^''' power-law dependence, which is to be expected if 
the toroidal current density obeys dl/dr oc r"^''*. 
The radial field roughly follows 



\Bn 



-y47rpo,diskc2 



0.4 I — 



(11) 



This radial dependence is close to the monopolar scaling of 
r~^, which is expected for a completely laminar flow where 
the mass inflow is conflned to a constant H/R. The monopo- 
lar radial fleld within the disk is associated with a. v — 
toroidal current within the core of the disk. Equations ([9]- 
llip show that the toroidal field dominates the other field 
components within the disk. 

This type o f solution is similar to the model proposed by 
iGammid (ll999D who described a solution such that within 
the plunging region the dimensionless radial magnetic fiux 



Fe6 = 



Fg, 



r^B"- 



V-27rr2poit'- \GMJ 



(12) 



is a constant function of radius, where we have temporarily 
reintroduced GM and c, B^ is in Gaussian units, and u"" is 
the radial 4-velocity. Here we report that this parameter is 
nearly constant throughout the entire accretion fiow out to 
r ~ 40rg with a value of 

Fe^ ^ 1.09, (13) 

which is the same as found bv lMcKinnev fc Gammid (|2004l ) 
for the plunging region in GRMHD simulations. The con- 
stancy of this parameter is consistent with a disk containing 
a radial field that is nearly monopolar per uint mass flux, 
as envisioned in the Gammie inflow model. 

This behavior of the accretion fiow is consistent 
with the fact that the radial dependence of within 
the GRMHD plunging region follows the Gammie in- 
flow solution |McKinnev fc Gammiell200^ . As described in 
iMcKinnev fc Gammie l|2004h . the thin disk Gammie solu- 
tion does well to model the GRMHD fiow apart from the 
lack of modelling the effects of pressure that lead to a non- 
zero radial velocity across the ISCO and the lack of a feature 
in the flow near th e ISCO. One primary difference is that 
the[G ammi 3 l|l999l ) model assumes Of is also constant along 
the radial field lines, while in the turbulent disk f^F ~ 
in the outer disk. 

In summary, despite the obvious turbulence in the bulk 
of the disk and the weak-disordered fleld in the corona, the 
time-averaged magnetic field averaged over the disk-|-corona 
has simpl e radial s caling s with a strong similarity to the 
model of iGammij (|l999l ). That model might be extend- 
able to apply to the entire accretion flow within a disk scale 
height. 



3.2 Comparison to the BP and ADAF Models 

The r''^^^ power-law scaling for the electromagnetic fleld is 
of particular interest because it occurs rather naturally in 
certain self-similar models in the literature. Here we check 
whether the assumptions made by these models applies to 
the GRMHD simulations. 

A popular model for disk winds is the 
iBlandford fc Pavnd (|l982l ) (BP) self-similar MHD wind 
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Figure 11. Radial profiles of the plasma rotational speed (v"^ , 
solid line), the magnetosonic speed (cma, dotted line), the sound 
speed (cg, short-dashed line), the Alfven speed [va, long-dashed 
line), and the radial 4- velocity (u'", dot - short dash line) for a 
GRMHD accretion disk simulation with a/M = 0.9375. All speeds 
are plotted in units of the Keplerian speed (i^k)- The angular and 
Keplerian speeds are in Boyer-Lindquist coordinates, while the 
others are given in comoving coordinates. The vertical lines indi- 
cate the radii of the ISCO and MBO. Note that the equipartition 
assumption, viz., Va ^ Cs ^ vk, is strongly violated. The solution 
at r > 20rg is still dependent on initial conditions and is excluded 
from the analysis. 



model in which they assumed that the Alfven speed at 
the equatorial plane scaled as the local Keplerian speed. 
Coupling this with the additional assumption that the 
density scales as p ~ r"^''^, they found that the magnetic 
field should scale as |&| oc r"^''*, which is consistent with the 
GRMHD flow as given by equation (O . The BP model also 
assumes that the disk is threaded by a large-scale organized 
field. Clearly this latter ass umption is b r oken, a s shown in 
figure [1] and as discussed in iHirose et al.l l|2004l ): lMcKinnevl 
|20"05a . 2006c,). 

More recently, the same scaling was obtained also 
in advection-domina t ed acc retion flow (ADAF) models 
iNaravan fc Yil 1 1994 [ 1995al ). Under the assumption of 
self-similarity, ADAF models naturally give p ~ ^.-3/2 
and pressure p ~ r~^^^ . Assuming equipartition between 

:as and magnetic pr essure, one again flnds |6| ~ ^-5/4 

Naravan fc Yilll995bl ). 

Figure [TT] shows the plasma rotational speed, magne- 
tosonic speed, sound speed, Alfven speed, and in-fall ra- 
dial 4-velocity for the fiducial GRMHD numerical model. 
The general relativistic generalization of the Alfven speed is 
given by 



(14) 



and the Keplerian speed is given by 



VK 



(15) 



lOrj, which 



Let us focus on the region interior to r 
has reached a quasi-steady-state (the region farther out 
is still sensitive to the initial conditions). Clearly, as Fig- 
ure [TT] shows, Va is not simply related to either vk or Cs . 
Thus, neither the BP nor ADAF assumptions are satis- 
fied. This is perhaps not surprising since self-similar mod- 
els assume a Newtonian gravity and so only apply far 
from the horizon. Figure [TT] shows that the magnetic field 
dominates over t he matter near the horizon, consistent 
with the results in De Vi llicrs. Ha wlev. fc KrolikI (|2003l ) and 
iMcKinnev fc Gammie ( 2004 ). This is also consistent with 
the fact that within the disk the ingoing fast magnetosonic 
surface is located at r ~ l.Srg, located between the ISCO 
and the marginally bound orbit (MBO). 

One important feature of black hole space-times is the 
ISCO, which for the simplest viscous thin disks demarcates 
where the fiuid plunges into the black hole. For our simu- 
lations of moderately thick disks with iif/i? ~ 0.26, we find 
that near the black hole within r < 3rg the fiuid plunges 
into the black hole with a power-law radial 4-velocity, 



-0.3 



(16) 



where is written in Boyer-Lindquist coordinates. For the 
fiducial model with H/R ^ 0.26 and a/M = 0.9375, this 
power-law plunging starts at approximately 



R 



plunge ' 



-RiSCO 1 + T7 



(17) 



This scaling is also consistent with the other H/R ~ 0.05 
model we studied and of course is trivially consistent with 
thin disk theory for which a thin Keplerian disk has the 
ISCO located at Risco = 2.044rg for a/M = 0.9375. Fur- 
ther studies can determine how general this fit is. However, 
notice that the fiuid begins to plunge toward the black hole 
at around r ~ 6rg, significantly further out than the ISCO. 
Thus there appears to be a factor of 2 ambiguity in iden- 
tifying the location where material plunges into the black 
hole. This GRMHD result has some not- well-defined bear- 
ing on recent measurements of black hole spin that used 
the o bserved spectra to estimate -Risco and so estimate 
a/M jShafee et al.ll2006l : [McClintock et al]|2006h . The in- 
going slow magnetosonic and Alfven surfaces are at r ~ 5r-g 
just inside where the fluid begins to plunge into the black 
hole. 

Within a scale-height and for all radii throughout the 
disk the mass accretion rate obeys 



M[r] « Const, 



(18) 



which indicates that for this particular model with a/M — 
0.9375 that mass-loss is not significant to the properties of 
the bulk flow. For more rapidly rotating black holes th e 
mass-loss can become signiflcant ([Hawlev fc Kroli^ |2006| ). 
but the radial dependence in the disk appears non-trivial. 
Since M[r] cc povT , then one expects a roughly constant 
proper density and indeed 



Po 



PO.disk 



r+ 



(19) 
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instead of po oc r^'^l"^ as one would expect in the BP or 
ADAF model. Notice that this defines po.disk- 

The growth in the Alfven speed is consistent with the 
fact that within r < lOrg the rest-mass density is nearly 
constant and the gas pressure is small and quickly diminishes 
at larger radii, following 

^L^O.Olf— ) (20) 

which implies that the enthalpy is small compared to the 
rest-mass density. An equally good fit has p oc r~^'^ . The 
comoving field energy is small compared to the rest-mass 
density, so that 

(va) OC OC r-^''^ (21) 
for r < lOrj. Interestingly, the de-correlated average 

^ '^\^\) (22) 

follows this dependence even more strictly than the direct 
space-time average of Va- This shows that space-time corre- 
lations are mild between the various sources of pressure. 

We thus conclude that it is purely an accident that the 
field and the current in the GRMHD solutions scale exactly 
as in the BP and ADAF models. Clearly the fact that Va oc 
r"^''* and vk oc r"^''^ means that ~ vk is not held and 
so the BP/ ADAF assumptions are violated in this region 
close to the black hole. This is found to be true for many 
models of the disk and a large range of black hole spins. It 
remains an open question as to what mechanism enforces 
the = 3/4 toroidal current density associated with the 
polar field. 
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Figure 12. Radial dependence of the effective viscosity param- 
eter Qmag as defined in eq. (19) for a GRMHD accretion disk 
simulation with a/M = 0.9375. The solid line corresponds to the 
choice P = 6'^/2, i.e., magnetic pressure alone, and the dotted line 
corresponds to P = /2 + pg, i.e., the total pressure. The dashed 
line shows the uncorrelated time-average of —2{V){b^)/{b'^). All 
quantities have been space-time averaged over a period of time 
covering several turbulent eddy time scales. The vertical lines in- 
dicate the radii of the ISCO and MBO. The curves show that the 
magnetic stress is significant inside the ISCO and that correla- 
tions among the field components are significant. 



3.3 Magnetic a viscosity parameter 

Accretion flows without magnetic fields are often assumed 
to be free of dissipation and torques within the ISCO. How- 
ever, it has long been understood that magnetic fields can 
drastically violate this assumption through the action of ex- 
tended fields that generate torques across the ISCO even 
without turbulent dissipation. While equation (|17p demon- 
strates that the gas pressure scale-height is related the ef- 
fective location of the ISCO, magnetic stresses may play 
some role in setting this scaling. Magnetic stresses may also 
play independent roles not at all modelled by a scale-height 
-based argument. 

An important parameter in understanding the pres- 
ence of stresses within the ISCO is the magnitude and 
radial scaling of the a viscosity parameter. Early esti- 
mates of the BZ power output depended upon the argu- 
ment that a was as determined in lo cal shearing box cal- 
culations (|Ghosh fc Abramowicj[l997l ). They also assumed 
that the field threading the hole was determined by a sub- 
equip artition between the field and gas pressure in the 
disk (| Ghosh fc Abramowicj|l997l : iLivio. Ogilvie. fc Pringld 
1 19991 ). Finally, they assumed that the electromagnetic power 
from the disk is not converted into material and ther- 
mal forms. Based upon these assumptions, they concluded 
that BZ power output was negligible and certainly weaker 
than the electromagnetic power of the disk. Here we check 
whether a behaves the same near the black hole as in local 
non-relativistic simulations. 



Local shearing box simulations of a small section of the 
accretion flow have suggested that Va ~ 2y^Cs, where a ~ 
0.01 — 0.1 is the usual dimensi onless viscosity in standard 
accretion disk models (jHawlev et al.. 1995). Global pseudo- 
Newtonian si mulations have shown th at a rises sharply in- 
side the ISCO (|Hawlev fc Krolikll200lh . In our GRMHD sim- 
ulations, the radial transport of angular momentum can be 
investigated by measuring the effective magnetic a 

amag ~ ( p ) , (23) 

where </!>'' — {0, 0, 0, 1} is the </> Killing vector associated with 
the axisymmetry of the system, is the comoving radial 
held strength, and 6^ is the covariant comoving 4-fleld. This 
form of Qmag is independent of the coordinate system for ax- 
isymmetric space-times. We choose P to be either the total 
pressure or just the magnetic pressur e b^/2. From the high- 
resolu tion fiducial model studied in iMcKinnev fc Gammiel 
(|2004l ). we compute the radial dependence of amag, inte- 
grated over a disk scale height {H/R ~ 0.26) and over the 
turbulent period of accretion. The angular integration is con- 
flned to the disk, so that the result is not influenced by the 
corona above the disk. The result is shown in Figure [I2l We 
see that Omag rises toward the horiz on, which is consisten t 
with the non-relativistic results of iHawlev fc KrolikI (|200ll ). 
The large Omag is associate d with a flux of angular mom en- 
tum from inside the ISCO (jMcKinnev fc Gammij|2004l '). 
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The plot also shows the uncorrelated time-average of 



^mag,uncorrclatcd 



(672) J 



(24) 



where the brackets denote a time-average. This quantity is 
nearly constant within the entire flow, showing that tempo- 
ral correlations between the quantities are significant. 

Since the value of a rises near the black hole, one might 
expect that the increased stress would be associated with an 
increased energy flux from the black hole and an associated 
increased luminosity from the disk near the ISCO compare d 
to standard thin disk theory l|Krolikl [l999l : lGammielll999l ). 
However, for these GRMHD models, which correspond to ra- 
diatively inefficient accretion flows (RIAFs), the energy per 
baryon accreted is similar to that for a thin disk with a fixed 
a an d no torques inside the ISCO l|McKinnev fc Gammid 
|2004| ). 

The fact that a ~ 1 near the black h ol e vio - 
lates the assumptions of Ghosh fc Abramowic j |l993) ; 
iLivio. Qgilvie. fc Pringlj l|l999^ . Also, black hole spin plays 
a nontrivial role in en hancing the magn etic field threading 
the black hole horizon (lMcKinnevll2005al ) and as discussed in 
section the field and gas pressure can be near equiparti- 
tion near the black hole. Finally, as discussed in section !?!^ a 
significant amount of the disk's electromagnetic power out- 
put is lost to material and thermal power. Thus, the BZ 
mechanism may still be responsible for the most of the elec- 
tromagnetic power output from accretion systems and may 
still account for the most powerful radio sources. These ef- 
fects should be considered when comparing the BZ and disk 
electromagnetic powers. 



4 LIMITATIONS 

The primary limitation of the present study is that the nu- 
merical models are axisymmetric. A 3D model may show 
that V = 3/4 is not generally chosen by the system. Com- 
parisons between 2D and 3D GRMHD simulations have 
shown reasonable qualitative and quantitative consistency 
toe ViUiers. Hawlev. fc KrolikI [iool iMcKinnev fc Gammiel 
l2004h . In particular, both show that the fiow partitions into 
a disk, corona, disk wind, and magnetized funnel region. 
Both give similar accretion rates of energy and angular mo- 
mentum per unit baryon. The primarily problem with ax- 
isymmetric simulations is that turbulence decays after some 
length of time. We avoid this problem by only making mea- 
surements during the turbulent period of accretion, so our 
results are unlikely to significantly change in 3D simulations. 

In the regime where radiation determines the energy 
balance in the disk, and so determines the disk thickness, ra- 
diative cooling should make little difference to these results 
since models with disk thicknesses between H/R ~ 0.05 to 
H/R ~ 0.3 show the same behavior. However, radiative ef- 
fects could also be dynamicall y important, such as throug h 
the photon bubble instability (|Aronslll992l : lGammielll998h . 
and it is unknown how this would change our results. 

The particular field geometry accreted can significantly 
change the mass-loading of the Poynting-dominated jet. Ac- 
creting a more random field leads to a larger coronal region 
that can extend all the way to the poles around the black 
hole, and no Poynting-dominated jet would form unless the 



black hole spin were sufficiently large that the toroidal mag- 
netic pressure exceeded the ram pressure of the coronal ma- 
terial. Preliminary models of Keplerian disks, with cooling 
to keep the thickness fixed, suggest that the accretion of an 
irregular field leads to no Poynting-jet formation for at least 
a/A/ < 0.94 (McKinney fc Gammie, in prep.). Despite the 
absence of the Poynting-jet, the toroidal current within the 
disk is still well-modelled by the = 3/4 solution. In the tan- 
gled field case, the corona simply fills the region that would 
otherwise have been occupied by the Poynting-dominated 
jet. 

All the models studied have disks that only extend out 
to r ~ 40r5 with a similar circularization radius. The re- 
sults of this paper, such as the radial scalings, only directly 
apply within r < lOr^ where general relativistic effects of 
the space-time play an important role. The results found 
here may not hold far from the black hole. For example, 
Qmag ~ Const. ~ 0.1 is expected far from the black hole. 
Also, while a more extended disk is not expected to affect 
the flow properties near the black hole, an extended disk 
may affect the Poynting-dominated jet at large radii. Future 
studies should investigate a more extended disk to model 
systems that have a large circularization radius. 



5 DISCUSSION AND CONCLUSIONS 

We have shown in this paper that GRMHD numerical 
simulations of magnetized accretion flows lead to a sim- 
ple angular-integrated toroidal current density of the form 
dl^/dr oc r~^/*, even though the accretion disk, corona, 
and outflowing wind are highly turbulent and chaotic. We 
showed that the poloidal fleld distribution in the jet is con- 
sistent with the force-free fleld solution for an r"^''* current 
distribution in a non-rotating equatorial disk. The fact that 
we obtained good agreement even with a highly simplified 
non-rotating model suggests that self-collimation by hoop- 
stresses may not be required within the jet. If hoop-stresses 
are ineffective, then force-balance requires the strongly mag- 
netized jet to be confined by the weakly magnetized corona 
and wind (see lLvnden-BellH2006l for a discussion of pressure 
confinement). By measuring forces across the corona-funnel 
interface, we found that force balance primarily requires the 
coronal pressure indicating that the corona is required to 
confine/coUimate the polar jet. Future work should study 
by what mechanism the confining coronal pressure generates 
a force-free field consistent with the r~^^^ toroidal current 
density. 

Because of the strong turbulence, at large radii the 
plasma and the magnetic field in the disk are locked to- 
gether so that both rotate around the black hole at roughly 
the local Keplerian angular frequency (Oif). For r < Svg, the 
behavior of the plasma and field qualitatively changes and 
depends strongly on whether the black hole spin is larger or 
smaller than a/M ~ 0.4. Near the horizon, the field angular 
velocity asymptotes to Qf ~ f^£r/2 in the case of rapidly 
spinning holes with a/M > 0.4, where Qh is the angular 
frequency of the hole, and to Qf ~ for a/A/ < 0.4. 

Although the disk is turbulent, the average magnetic 
field strength in the disk varies smoothly as |6| oc r"^''*, and 
the individual components of the field (/?'", , and B"^) also 
have simple power-law scalings. The scaling of |&| is similar 
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to that assumed in the Blandford & Payne (1982) and ADAF 
(Narayan & Yi 1995b) models, but we showed that the simi- 
larity is purely accidental. Those models assume Newtonian 
gravity and that gas and magnetic pressures scale similarly 
with radius. Gravity near the black hole is obviously not 
Newtonian, and in the simulations the magnetic pressure 
rises more rapidly then the gas pressure toward the black 
hole. Thus, the physical reason for the particular scaling of 
the field with radius i s unclear . Ther e is some resemblance 
to the inflow model of lGamrtii3 (|l999') , but the construction 
of a simplified 1-D model is left for future work. 

The dominance of the magnetic field near the horizon 
is associated with an increased magnetic stress that implies 
a magnetic viscosity parameter a ~ 1, while a ~ 0.1 far 
from the black hole. The increased magnetic stress leads 
to enhanced angular momentum transport, but the energy 
per baryon accr eted remains similar to that from thin disk 
theory (see also iMcKinnev fc Gamrtii^ 120041 '). Prior studies 
that estimated the BZ power of black holes and the disk 
power output assumed an a b ased upon local shearing-box 
models that obtained q ~ 0.1 l|Ghosh fc Abramowicj|l997l : 
iLivio. Qgilvie. fc PringQl 19991 ). These studies also assumed 
that the field strength near the black hole would be set by 
sub-equipartition arguments, but we find that the field near 
the black hole is a non-trivial function of black hole spin 
(jMcKinnevI l2005al ) . They also assumed that the disk was 
threaded by an organized field, but we find that turbulence 
disorganizes any ordered field threading the disk. The lack 
of an organized field threading the disk is associated with a 
significant amount of electromagnetic power from the disk 
being converted into material and thermal forms. In light of 
all these facts, estimates of the BZ power (and its compari- 
son to the disk power) should be reconsidered. 

The basic picture of the disk-jet coupling that emerges 
is that turbulence driven by the MRI leads to toroidal cur- 
rents that in a force-free state would drive an organized 
poloidal field from the disk. However, turbulence dominates 
the accretion flow and large-scale modes lead to angular mo- 
mentum transport that continuously advects equatorial field 
through the disk into the hole and transports field at higher 
latitudes to larger radii. This keeps the corona free of or- 
dered fields. The disk-corona interface is a site for dissipa- 
tion of disk field through reconnection. The energy from the 
dissipation drives a hot coronal wind off the disk. The coro- 
nal wind is thus quite different from MHD models of disk 
winds (e.g., Blandford & Payne 1982), and is more akin to a 
thermally driven wind. The pressure of the corona confines 
the Poynting-dominated jet, since without the corona the 
field gradients at the funnel-corona interface could not be 
supported. 

The geometry of the accreted magnetic field plays a 
crucial role in controlling the production of the Poynting- 
dominated jet. Quite similar jets are obtained in GRMHD 
simulations that start from a variety of different initial 
configurations of the magnetic field: uniform vertical field, 
poloidal loop of magnetic field in the disk, multiple loops of 
alternating poloidal direction. Even though the latter two 
models have zero net vertical field, nevertheless, they end 
up with a net vertical fiux through the black hole and jet. 
The mechanism by whi ch this happens is describe d in the 
discussion of Model B in llgumenshchev et al.l (|2003l ) and see 
iNaravan et al] (|2003l ). However, for models initialized with 



a mostly disorganized field, the Poynting-dominated jet is 
weaker or absent because mass continuously loads the polar 
region. Despite the lack of a simple Poynting-dominated jet 
in such models, there is still a disk wind containing a disor- 
ganized field and the currents and field strengths within the 
disk still follow the same power-law dependencies. 

Although we obtain good agreement between the mag- 
netic field geometry in a simple force-free current sheet 
model and the poloidal magnetic field configuration in the 
Poynting-dominated jet in GRMHD models, we note that 
we have not solved for the force-free toroidal structure of 
the field in the jet. Since our force-free model has no rota- 
tion, B*^ — and so there is no Poynting flux. Also, the 
Lorentz factor cannot be estimated from our non-rotating 
force- free solution. To model these quantities we must con- 
sider force-free models in which the disk and field are rotat- 
ing (e.g., with a Keplerian profile). It would be interesting 
to see how well simple rotating force-free models can repro- 
duce the toroidal structure and acceleration of Poynting- 
dominated jets found in GRMHD simulations. This is the 
topic of a companion paper (McKinney & Narayan 2006). 
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APPENDIX A: EVOLUTION EQUATIONS 

The GRMHD equations of motion are used to study mag- 
netized accretion disks in the gravitational field of rapidly 
rotating blac k holes described by th e Kerr metric using the 
HARM code (Gammio ot al."2003a|}_ with improvements de- 
scribed in McKinney (2006c); No ble et all l|200d '). The Kerr 
metric is written in Kerr-Schild coordinates, such that the 
inner-radial computational boundary can be placed inside 
the horizon and so out of causal contact with the flow. The 
Kerr metric in Kerr-Schild coordinates and the Jacobian 
transformation to Boy e r-Lin dquist coordinates are given in 
iMcKinnev fc Gammid l|2004h . 

Boyer-Lindquist coordinates are not chosen because it is 
difficult to avoid interactions between the inner-radial com- 
putational boundary and the jet. The coordinate singular- 
ity at the event horizon in Boyer-Lindquist can be avoided 
by placing the inner-radial computational boundary outside 
the horizon. However, Poynting-dominated flows have waves 
that propagate outward even arbitrarily close to the event 
horizon. Using Boyer-Lindquist coordinates can lead to ex- 
cessive variability in the jet since the ingoing superfast tran- 
sition is not on the computational grid, and then the details 
of the boundary condition can signiflcantly impact the jet. 
Numerical models of viscous flows h ave historically had re- 
lated issues (see discussion in, e.g., IMcKinnev fc Gammid 
l2002ll . 
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Al GRMHD Equations of Motion 

The GRMHD notation follows iMisner et all (| 19731 ). here- 
after MTW. A single-component MHD approximation is as- 
sumed such that particle number is conserved, 

(Po^t'');^ = 0, (Al) 

where po is the rest-mass density and is the 4-velocity. 
A 4-velocity with a spatial drift is introduced that is unique 
by always being related to a physical observer for any space- 
time and has well-behaved spatially interpolated values, 
which is useful for numerical schemes. This 4-velocity is 

u' = u" — "frf , (A2) 

where 7 = —u'^rja- The additional term represents the spa- 
tial drift of the zero angular momentum (ZAMO) frame de- 
fined to have a 4-velocity of r;^ = {—a, 0,0,0}, where a = 
Bind so u' = 7/0:. One can show that 7 — (l+q^)^^^ 
with = QijU^u' . 

For a magnetized plasma, the energy-momentum con- 
servation equation is 

r'"';. = (r^A + rE^M);. = o. (A3) 

where T*^" is the stress-energy tensor, which can be split 
into a matter (MA) and electromagnetic (EM) part. In the 
fluid approximation 

^^A = (Po + ^^^)^V-fp,P^^ (A4) 

with a relativistic ideal gas pressure pg — (7 — l)ug, where 
Ug is the internal energy density and the projection tensor 
is P''" = g'^" + u'^u" , which projects any 4- vector into the 
comoving frame (i.e. P^^Uf^ = 0). 

In terms of the Faraday (or electromagnetic field) tensor 



and 



-'em 



af3, 



(A5) 



which is written in Heaviside-Lorentz units such that a fac- 
tor of 47r is absorbed into the definition of F''", where the 
Gaussian unit value of the magnetic field is obtained by mul- 
tiplying the Heaviside-Lorentz value by \/47r. The induction 
equation is given by the space components of -F = 0, 
where F^" = |e^'"^^-F'reA is the dual of the Faraday ten- 
sor (Maxwell tensor), and the time component gives the 
no-monopoles constraint. Here e is the Levi-Civita tensor, 
where e'^'^^^ — — [fiuXS] and [pi^AiJ] is the completely 
antisymmetric symbol. The comoving electric field is defined 
as 

1 



jj.h'kX 



Uv Fxk = Tjf, 



(A6) 



where r) corresponds to a scalar resistivity for a comoving 
current density j'^ — J^P"'^. The comoving magnetic field 
is defined as 



6-- = u^F^" = -e^-"^^ 



u^Fk\ 



(A7) 



The ideal MHD approximation, 77 = e*" = 0, is assumed, and 
so the invariant e'^b^ = 0. Since the Lorentz acceleration on 
a particle is = qe'^, then this implies that the Lorentz 
force vanishes on a particle in the ideal MHD approximation. 
Since e^u^ = fe'^it^ — 0, they each have only 3 independent 
components. One can show that 

b^u" ^b-'u", (A8) 



(A9) 



so that the electromagnetic part of the stress-energy tensor 
can be written as 



-'em 



-{u'^u'' + Pi"')-b''b" 



(AlO) 



The other Maxwell equations. 



(All) 

define the current density, J*^, but are not needed in the 
ideal MHD approximation for the evolution of the matter or 
the magnetic field. 

For numerical simplicity, another set of field vectors are 
introduced, such that = F and Ei = Fu/^—g. The 
two 4-vectors e'^ and b'^ and the 3-vectors B* and Ei are 
just different ways of writing the independent components 
of the Faraday or Maxwell tensors. Equation (|A7|I implies 
6* = B''Ui and fe' = (B' + u^b^)/u*. Then the no-monopoles 
constraint becomes 



{^B^),, = 0, 

and the magnetic induction equation becomes 

iV^B^),t = -iV^ib^u' -b'u^)),j 
= -iV^iBV -B'v^),, 



(A12) 



(A13) 



where = u^/t 



Ei = £i 



-eijkV^B'' = — V X B is the 



EMF, and e'-''° is the spatial permutation tensor. The above 
set of equations are those that are solved. A more complete 
di scussion of t he relativistic MHD equations can be found 
in lAnild (|l989l ). 

A2 Stationary, Axisymmetric Constraints 

We now write down the Faraday tensor in terms of a vector 
potential A^, where _F)i^ — Av^^ — Afj,^v If the field is axisym- 
metric (9^ — > 0) and stationary [dt — > 0), then evaluating 
the condition = one finds that 



A^^eAt^r - At,eA, 
It follows that one may write 

At,0 At,r 



0. 



A. 



<t>,B 



A(t 



= -Q.F{r,e) 



(A14) 



(A15) 



F'^ 



where Q.F(r,6) is an as-yet-unspecified function. It is usu- 
ally interpreted as the "rotation frequency" of the electro- 
mag netic field (this is Ferra ro's law of isorotation; see e.g. 
I Frank. King, fc Rainel |2002| . 59.7 in a nonrelativistic con- 
text). This ca n also be written as Of = Ftr /Frd = Ftg/Fe,p. 
As shown in iMcKinnev fc Gammid (|2004l ). one can then 
write Ffi^v in terms of the three free functions Q,f,A^, and 
B"^ , the toroidal magnetic field: 

(A16) 
(A17) 
(A18) 
(A19) 



Ftr = 


— Frt — ^pA^^r 


Fte = 


—Fet = ^.pA^fi 


Fre = 


-Fer = V^B'>' 


Fr4, ■- 
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Fe. 



(A20) 



with all other components zero. Written in this form, 
the electromagnetic field automatically satisfies Maxwell's 
source-fi-ee equations. Notice that A^^g = ^~gB^ and 



A 



APPENDIX B: REDUCING THE GRMHD 
SOLUTION TO A CURRENT SHEET 

Once we have trivialized the full GRMHD accretion disk 
into a toroidal current sheet, then we can obtain the effec- 
tive stationary, axisymmetric poloidal magnetic field from 
equation HA11|) . For an axisymmetric, stationary solution 
equation (|A11|) gives that 



(Bl) 



where Fe,/, 



-gB'', Frtj, = A^ 



pr<l= ^ gt''-g''"l>Ff,^. The equation in Boyer-Lindquist be- 



+ (\g\g'>'B^{Qpg'^-g^^)),e, 
or in terms of the vector potential and flp , 



(B2) 



+ 
+ 



{V^g'">A^A^pg'^-g^^)),o 



(B3) 



In either Boyer-Lindquist or Kerr-Schild coordinates, 
the linear partial differential equation (PDE) for a/AI — 
is 

' -2M ■ 



- sin 9 J''' = 



and for a current sheet, 



2M 



-A 



4>,e 



(B4) 



Aj, 



+ 



,A. 



4>,0 



(B5) 

This equation is an elliptic partial differential equation for 
r > 2M and is hyperbohc for r < 2M. The quantity A''* is 
the surface current density on the equatorial current sheet. 



Bl Vacuum Solutions in Schwarzschild space-time 



Equation (IB5|l can be solved to find the complementary so- 
lution by setting the quantity on the right to everywhere 
except on the current sheet at 9 — n/2 and assuming that 
the solution is separable such that 



A^{r,9) = R{r)e{9). 



(B6) 



Since the PDE is second order, there are in general two free 
functions. The radial ordinary differential equation (ODE) 
can be written in closed form in terms of generalized hyper- 
geometric function or solved numerically for AI ^ 0. In the 
limit of a; oo, one finds a complementary function of 

R{r) = Cor-' + Cir+\ (B7) 

where Co and Ci are arbitrary constants and / ~l/2. 
These are the standard spherical radial eigenfunctions. The 



angular complementary function can be written in terms of 
generalized hypergeometric functions or alternatively writ- 
ten in terms of the associated Legendre functions of the first 
(P) and second (Q) kind, where 



6(6') = Dol sin6i|P,^(cos6') -FDil sin6'|(3i^(cos( 



(B8) 



where Do and Di are different arbitrary constants and where 
I ^ —1/2 are linearly independent (i.e. P^;_i = Pt)- This 
form is a mixture of t wo hypergeometric functions f or each 
I (see, e.g., chpt. 8 in Abramowitz fc Stegun|[l973 and in 
iGradshtevn fc Rvzhiklll994 ). These functions are just one 
of the vector spherical harmonics. 



B2 Constraints 

This general solution is a sum of any coefficients for all I, but 
has no constraints to avoid divergences (e.g. monopoles or 
divergences in the physical field strength) on the coordinate 
singularities. The only constraint required is that F'^'^F^,^ 
remains finite, where in Boyer-Lindquist coordinates with 
a/M = and B'' = 0, 



F^''F^, = 2{{F^^Ygrrg^^ + {F 



(B9) 



and so to avoid divergences at the polar axes one requires 
A4, oc 6*^+" 4- Const, for n > 0. There is no requirement on 
the horizon and no constraint on divergence need be placed 
at r = since it is a physical singularity. 

The solution given by equation ([BSP with sin(S) multi- 
plied by the associated Legendre functions of the first kind 
gives A^ (X 9^ -\- 0{9^) for any Z, while the term with the 
second kind gives cx Const. + 0{9'^) for all I. Hence, 
the only valid solutions are of the first kind with any ra- 
dial solution, except for angular solutions combined with a 
radial function that does not depend on radius. The single 
nontrivial example of this exception is the monopole solu- 
tion with I = —1 and Co = Do = giving A^ — — cos 9. 
The remaining solutions require the first kind with Di — 0. 
The monopole solution is only a result of the limit to the 
open set around I = —1 for Ci 7^ or Z = for Co 7^ 
for the associated Legendre function of the first kind. This 
issue regarding the monopole solution can be considered a 
pathology of using the Legendre functions that is not man- 
ifested in the hypergeometric form of the solution, where 
/ — naturally generates the monopole and paraboloidal 
type solutions. 

B3 Currents 



The current density given by equation (lAllfl can be inte- 
grated to determine the net toroidal current enclosed within 
the volume between ro and r, as given by equation [5] For a 
model with a current sheet located at 6* = 7r/2 with 

J-f" ^ K'^5{9-7v/2), (BIO) 

where is the surface current density, then 

^ = V^ie = 7v/2]Kt (Bll) 
For a general power-law toroidal current per unit radius 

with 

dl^ _ _C_ _ _C_ 



(B12) 
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the enclosed current (for n 7^ 0) is 



(B13) 



where r2 denotes some radius for which ro < r < r2 and 
for which the enclosed toroidal current I<i,(r2) is known. If 
a/M = and in spherical polar coordinates, then \/^~g[6 = 
it/ 2] = r^, so that if the toroidal surface current density is 

C C 



(B14) 



B4 Split-Monopole Field 



The split-monopole solution for a/M = is a solution where 
the current sheet at 6 = n/2 has 



J^ = ^5(^-^/2), 

so that 

dl,p _ C 
dr ' 

which gives a total enclosed current of 

T - ^ 

^ J 1 ' 

ro r 

where ro < r are the inner- and outer- radial positions. The 
split-monopole has n = 1 and z/ = 0. 

The linearly independent particular solution for the 
split-monopole vector potential (for any M) is 



(B15) 



(B16) 



(B17) 



A 



(split) _ j -C COS e < n/2 

+Ccose(l>f, 9>7v/2, 



(B18) 



where (p'^ = {0,0,0,1} is the (j> Killing vector. The above 
gives that = B"* = and that 



B" = 



+ ^ e<n/2 
e>-n/2, 



(B19) 



This solution has field lines that have an opening angle 
following 

Bj oc r°. (B20) 



B5 BZ Paraboloidal Field 

The paraboloidal field solution for a/M = is a solution 
where the current sheet &t 9 = n /2 has 



so that 



dl^ _ C 
dr r 



(B21) 



(B22) 



Notice that the total enclosed toroidal current would diverge 
for a disk of infinite radial extent, i.e., 



/^ = Clog(- 

\ro 



(B23) 



where ro < r are the inner- and outer- radial positions. The 
paraboloidal solution has n = and v = 1. 



The linearly independent particular solution for the BZ 
paraboloidal field vector potential is 



^(para) ^ ] +3 (r, 6^) <j!>^ 6 < 71/2 

+3(r,7r-e)0^ 9>n/2, 



(B24) 



where g(r,0) = %[rf. + 2M/+(1 - ln/+)], f+ = l + cos9, 
/_ = 1 — cos9, and M is the mass of the black hole. Thus, 
B'f' = and 



B" = 



+ 



C(r+2Ml-nf_) 



and 



2r2 



— . Ccot{e/2) 

2P 



>^/2, 



> V2, 



(B25) 



(B26) 



This solution has field lines that have an opening angle 
that at large radii approximately obeys 

9j oc r"° ^ (B27) 



B6 Constructing Current Sheets by Splicing 
Source-Free Solutions 

While in general it is difficult to find for arbitrary source 
functions, one can construct source functions corresponding 
to equatorial current sheets by splicing a single equatori- 
ally asymmetric vector potentials G{r,9). Then the general 
solution with a current sheet is 

A^{9) = G{9){l-H{9-Tv/2)) + G{n-9)H{9-n/2), (B28) 

for either G in the range of ^5 ^ ^ n/2 or n/2 !^ 9 ^ n, 
where H is the Heaviside function. 

First the simple a/M = M — equations are consid- 
ered. To construct a radial power-law current density of the 
form given by equation HB14|I so that the source function 



the radial complementary or particular 



functions must satisfy 



R{r) oc r" 



(B29) 



at large radii (r ^ M) as demonstrated by plugging this 
form of Atf, into equation (|B4p . This means that one must 
choose either 



with Co 7^ and Ci 



I = -V 
or one must choose 
l = v-l 



(B30) 



(B31) 



with Ci 7^ and Co = 0, where / ^ —1/2. Hence, for v > 
1/2, only the second choice leads to a real solution. For < 
1/2, only the first choice leads to a real solution. For f = 1/2 
both give / = —1/2. Thus, for fixed i/, there is a unique I 
that gives a single allowed R{r) oc r" and 0{9). 

For example, for v — 1 one has that 1 = 0, and then the 
most general solution for M = is 



A4, = (co + cir)(di -I- d2 coat 



(B32) 



which after forcing A^ oc 9"^ + Const, for ^ 1 near 9 = 0, 
one has that 



which for co 



^0 = (co -hcir)(cos6'- 1), (B33) 
gives paraboloidal solution in the upper 
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hemisphere, ci — gives the monopole solution, while com- 
binations give a mixture of monopole and paraboloidal solu- 
tions. Another example is Co = Di — and 1 — 1/2 giving 
a decollimating field geometry. 

The GRMHD numerical solutions are associated with 
models with v — 3/4, for which one must choose I = —1/4 
and Co = 0. For the solution to satisfy regularity on the axis 
near 9 — 0, one must set 

Do / and Di = 0. (B34) 

This solution is quite similar to the paraboloidal solution, 
but slightly less coUimated, as expected. This solution has 
field lines that have an opening angle that approximately 
follows 

6'i oc r-°■^^^ (B35) 

This force-free mo del is also used in 
iMcKinnev fc NaravanI (|2006[ l to find force-free solutions 
with arbitrary M and a/ M using a general relativistic force- 
free e lectrodynamic code (|McKinnevll2006bl : iNaravan et al.l 
l2006l ). 
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